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ABSTRACT 

We study the compact binary population in star clusters, focusing on binaries 
containing black holes, using a self-consistent Monte Carlo treatment of dynamics and 
full stellar evolution. We find that the black holes experience strong mass segregation 
and become centrally concentrated. In the core the black holes interact strongly with 
each other and black hole-black hole binaries are formed very efficiently The strong in- 
teractions, however, also destroy or eject the black hole-black hole binaries. We find no 
black hole-black hole mergers within our simulations but produce many hard escapers 
that will merge in the galactic field within a Hubble time. We also find several highly 
eccentric black hole-black hole binaries that are potential LISA sources, suggesting 
that star clusters are interesting targets for space-based detectors. We conclude that 
star clusters must be taken into account when predicting compact binary population 
statistics. 
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1 INTRODUCTION 

The inspirals and mergers of compact binaries where both 
members are neutron stars (NS) or black holes (BH) are 
some of the most promising sources for the current and next 
generation of ground based gravitational wave (GW) detec- 
tors (LIGO, Virgo, GEO 600, TAMA 300). NS-NS binaries 
are expected to be the most plentiful merger species in the 
frequ e ncy regime of ground b ased detectors (jAbbott et al.l 
120051 : IBelczvnski etafl l2007h . however BH-BH binaries 
are more massive and thus can be detected at larger 
distances, as far as the Virgo cluster ( Abbott et alj [20061 ) 
for the current generation of gravitational wave detectors 
and up to cosmological distances for the next generation. 
Compact binaries with longer periods and including white 
dwarfs (WD) may also be detectable in the low-frequency 
(5 x 10~ 5 — 1 Hz) LISA (Laser Interferometer Spa ce 



Antenna) frequency band dHils. Bender fc Webbink 1990l: 
[Benacquista 2001; Nclc mans. Yungelson fc Portegies Zward 
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120011 : IBelczvnski. Benacquista fc BulikI |201Ch . WD- WD 
binaries will be the most plentiful stellar mass LISA 
sourc es and are expected to produce confusion limited noise 
(e.g. lEvans. Iben fc Smarn Il987l: iHils. Bender fc Webbinkl 
1990l: Nelemans. Yungelson fc Portegies Zwar ij l200ll : 
Timpano. Rubbo fc Cornish! 1200(1 iRuiter et ail 120071 ) 
whereas the less common NS-NS, NS-BH, and BH-BH 
binaries are potentially resolvable. 

In order to make predictions for GW event rates, the 
population of compact binaries in the universe must be un- 
derstood. While the population of NSs in the local uni- 
verse can be constra i ned by observ ations of pulsars (e.g. 
iKalogera et al. 1 120011 : lLorimedl2005h , BHs cannot be ob- 
served directly and their properties can only be constrained 
by modelling. There have been several studies carried out on 
the compact binary population in the galactic field where 
stellar and binary evolution proceeds in isolation and can 
be modelled using s i mple population synthesis. In particu- 
lar [BeIczxnskret^L| l|2007h predict a detection rate for the 
advanced LIGO detector of ~ 20 NS-NS mergers yr _1 , only 
~ 2 BH-BH mergers yr _1 , and ~ 1 NS-BH merger yr _1 . 



2 J. M. B. Downing, M. J. Benacquista, M. Giersz, & R. Spurzem 



Thus the detection rate in the galactic field should be dom- 
inated by NS-NS mergers. iBelczvnski. Benacquista &: Bulikl 
(2010) have performed similar calculations for the galactic 
field in the LISA band. Depending upon the assumptions 
made about the probability of mergers during common en- 
velope evolution they find 2 — 6 resolvable NS-NS binaries 
and — 5 resolvable BH-BH binaries. This implies that, al- 
though rare, both NS-NS and BH-BH binaries may appear 
as resolved stellar mass sources in the LISA band. Overall 
NS-NS binaries will dominate the field detection rate for 
ground-based detectors while a detection in the LISA band 
is possible but unlikely. 

In star clusters the situation is rather different. Here 
interactions between stars and binaries are common (e.g. 
lHeggieiri975h and can affect the final outcome of binary evo- 
lution. Such interactions can form new binaries from single 
stars, exchange binary members and field stars and can re- 
duce or increase the period of existing binaries. In order 
to maintain energy equipartition, interactions between par- 
ticles of different masses tend to accelerate the lowest-mass 
particles to the highest velocities (ISpitzerlll987T ) . As a conse- 
quence low-mass objects are the most likely to escape during 
few-body encounters. Therefore few-body interactions tend 
to introduce massive objects into binaries. Massive objects 
also tend to sink to the center of star clusters where the stel- 
lar density is highest ( mass se gregation, also a consequence 
of energy equipartion, ISpitzeil Il987l ) and thus massive ob- 
jects are the most likely to experience dynamical interac- 
tions. Since BHs rapidly become the most massive objects 
in star clusters due to stellar evolution they will be partic- 
ularly strongly affected by dense stellar environments and 
are very likely to be exchanged into binaries. Therefore star 
clusters are predicted to form BH-B H binaries rather effi- 
ciently (jSigurdsson fc Phinnevlll993l ) and may significantly 
enhance the BH-BH merger rate in the universe. 

Several authors have inv estigated this possibility us- 
ing v arious approximations. iGiiltekin, Miller fc Hamilton! 
l|2004 ) and lO'Learv etail (|2006l ) have simulated the for- 
mation of intermediate-mass black holes (IMBHs) assuming 
that the BHs in the cluster are completely mass-segregated 
and interact only with each other. In this situation, where 
the BHs interact very strongly, the authors resolve these 
encounters with explicit few-body integration. The interac- 
tions can lead to the efficient formation of BH-BH binaries 
but also tend to destroy or eject those already formed. These 
authors conclude that gravitational wave mergers can oc- 
cur within young star clusters but BH-BH binaries will be 
destroyed rather efficient ly and the populatio n will be de- 
pleted within a few Gyrs. lO'Learv et al.l (|200q ) in particular 
calculate a star cluster BH-BH detection rate for advanced 
LIGO of 1 — 10 merges yr" 1 , up to 70% of which actually 
occur in ejected BH-BH binaries and thus take place in the 
galactic field. Neither set of autho rs include stell a r evo lution 
in their simulat i ons. B y contrast llvanova et all (|2008T l and 
ISadowski et al.l (|2008T ) have conducted studies of the com- 
pact binary population in star clusters using stellar evolu- 
tion prescriptions and simplified two- zone models of cluster 
dynamics. The two-zone models assume that the BHs and 
BH-BH binaries remain in dynamical equilibrium with the 
rest of the cluster and do not strongly mass-segregate. Thus 
the density of BHs and BH-BH binaries remains low and 
they do not interact with each other nearly as frequently as 



in the previous models. ISadowski et al.l (|2008T l in particular 
find no NS-N S mergers but a much higher rate of BH-BH 
mergers than lO'Learv et alj (|2006l ). They calculate a de- 
tection rate of 25 — 3000 mergers yr _1 for advanced LIGO 
even though thei r treatment of the fe w-body interactions 
is the same as for lO'Learv et al. (2006). This is because al- 
though there are fewer interactions that create BH-BH bina- 
ries there are also fewer interactions that destroy them. This 
highlights the importance assumptions about global cluster 
dynamics can have on detection rates. The only study with 
fu ll treatment of both dynamics and s tellar evolution is that 
of lPortegies Zwart fc McMillan! (2000) who conducted small 
(iV ~ a few 10 3 ) direct N-body simulations and showed that 
BH-BH binaries are quickly ejected from star clusters due 
to strong few-body inte r action s. This seems to confirm the 
model of lO'Learv et al] (|2006l ) but the simulations are too 
small for any strong conclusions to be drawn. It seems that 
star clusters can significantly enhance the rate of BH-BH 
mergers in the universe however by exactly how much de- 
pends strongly on the dynamical assumptions made. 

All of these simulations resolve the few-body interac- 
tions but either use very simplified models for the global 
cluster dynamics or have values of N too small to repre- 
sent GCs. In this paper we use a Monte Carlo code to self 
consistently model the dynamics of globular clusters over a 
range of metallicities, binary fractions, and initial concentra- 
tions. We hope to constrain which dynamical assumptions 
are most likely to produce accurate results for gravitational 
wave event rate predictions. We will focus only on the com- 
pact binaries that can be found within the cluster during 
its evolution, leaving a detailed discussion of the population 
of escapers to a future paper. In §[2] we briefly describe the 
Monte Carlo code and some of its features and limitations. 
In §0we describe our initial models. In §fj]we present the 
results of our simulations. In §[5] we present predictions for 
LISA detections. We discuss our results in §[6] and conclude 
in § El 



2 NUMERICAL METHODS 

Monte Carlo star cluster simulations use Monte Carlo in- 
tegration of the theory of two-body relaxation in order to 
approximate the evolution of GCs. Assuming a spherically 
symmetric potential, the orbit of each centre of mass (sin- 
gle star or binary) in the cluster at any instant can be de- 
fined by its energy, E, and angular momentum vector, J. 
Changes in E and J due to the surrounding stars can then 
be calculated by an appropriate choice of random scattering 
events from the theory of two-body relaxation. In this way 
the dynamical evolution of the star cluster can be simulated 
self-consistently. A position for each star can be defined us- 
ing a time-weighted average over the orbit defined by E 
and J and then the probability of encounters between stars 
can be calculated. Unlike other approximate methods for 
calculating star cluster evolution, each centre of mass is ex- 
plicitly included in the simulation. Therefore it is relatively 
straightforward to incorporate special prescriptions for in- 
dividual astrophysical events such as few-body interactions 
(which are not covered by the Monte Carlo approximation) 
and stellar evolution. 
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2.1 The Monte Carlo Code 



We use a Henon-t ype code rtHenorj 197ll ) inc orporating the 
impro vements of IstodolkiewiczT l 1982) and IStod61kiewic3 
ll 19861) for both global and binary dynamics as described 
bv lGiers3 (|l998l ). The code includes prescriptions for three- 
body binary formation and for binary-single and binary- 
binary encounters. The probability for three-body binary 
formation interactions are calculated between all adjacent 
st ars at each ti mestep according to equations (7) and (8) 
m iGiersd (|200ll ). New energies, vel ocities, and o rbital pa- 
rameters are calculated according to lGiersa l|l998l) . Binary- 
single and binary-binary i nteraction probabili t ies ar e cal- 
culated in a similar way to iGiersz fe Spurzeml (|2003l ). The 
outcome of the binary-bin a ry in teractions follow the pre- 
scriptions of IStodolkiewiczl (Il986l) which are in turn based 
on the numerical experiments of iMikkoial (|l984l ). Exchange 
interactions, where one binary member can be exchanged 
for a field star or the member of another binary, are al- 
lowed during binary-single and binary-binary interactions. 
The most common exchange outcome is a light star be- 
ing exchanged for a massive star. The probability of an ex- 
change interaction for each star is given by equation (17) of 
IHeggie. Hut fc McMillanl |l996|). Tidal t runcation is treated 
in a simplified way using the approach of Baumgardtl (120011 ) 
(|Gierszll200ll : IGiersz. Heggie fc Hurlevll2008h . 

Single and binary stellar evolution are simulated us- 
ing the Single Stellar Evoluti on (SSE) and Binar y Stel - 
lar Evolution (BSE) recipes of IHurlev. Pols fc Tout! (|2000|) 
and IHurlev, Tout fc Pols! |2002i ) jGiersz. Heggie fe Hurlevl 
2008). These recipes include a full treatment of binary and 
stellar evolution from the zero age main sequence to the de- 
generate remnant for a variety of stellar masses and metallic- 
ities. Of particular interest to our work are the natal velocity 
kicks applied to NSs and BHs due to asymmetric super- 
nova explosions (|Lvne fc Lorimerlll994l ). Velocity kicks are 
applied to all NSs and BHs at birth and are drawn from a 
Maxwellian velo city distribution with a dis persion of ~ 190 
km s _1 based on lHansen fc Phinnevl (|l997i )'s proper motion 
samples of NSs. The kick velocity of the BHs is then reduced 
in proportion to the mass of accr e ted m aterial as described 
in iBelczvnski. Kalogera fc Bulikl (|2002l ). The survivability 
of the binary is calculated as described in IBelczvnski et all 
(2006). A simplified grav itational wav e inspiral timescale in 
the weak field limit (e.g. |Peterslll964l ) is also included. 

The code has been compared to direct N-body 
simulations and produces excellent agreement b etween 
globa l dynamical properties for both s ingle ( Gierszl 



19981) and multi-mass c ases jGierszl l200ll . IGiersz! 



2006). 



Giersz. Heggie fc Hurlevl (|2008l ) have shown the code com- 



pares very well with direct N-body simulations when stellar 
evolution is included. The code is also able to re-produce 
several observed physical parameters, such as the surface 
brightness profiles an d luminosity functions, of t he ob served 



star clusters M67 (IGiersz. Heggie fc Hurlevl 120081). M4 
(IHeggie fc Gierszl l2008h and NGC 6397 ([Giersz fc Heggid 
2009). 

An attractive feature of Monte Carlo simulations is that 
the computation scales with 0(N 1 ) — 0(N 2 ) rather than the 
0(TV 3 ) of direct TV-body codes. Furthermore including bina- 
ries does not greatly decrease the performance. This means 
that a simulation of TV m 5 x 10 5 — 10 6 bodies with 50% bina- 



ries can be carried out on a single fast processor in the order 
of hours to days rather than the weeks to months required for 
direct TV-body simulations. At the same time, unlike other 
approximate methods of calculating star cluster evolution, 
the Monte Carlo simulation manages to provide the same 
star-by-star information produced by direct TV-body codes. 
Thus the Monte Carlo code can be used for large parameter 
studies in a short space of time where the details of indi- 
vidual stellar events, such as inspirals and mergers, are of 
interest. 



2.2 Limitations 

The approximations used in the Monte Carlo code introduce 
limitations that may affect our results. In particular the BHs 
may be sufficiently massive compared to the rest of the sys- 
tem that the become "Spitzer unstable" ((Spitzcr 1987j) and 
form a decoupled small TV subsystem in the cluster core. This 
subsystem interacts only with itself and in such a situation if 
TV becomes small enough the distinction between large and 
small angle scattering breaks down. The Monte Carlo ap- 
proximation relies on the fact that these scales may be sep- 
arated (the cluster can be divided into a near and far zone) 
and may becom e unre liable for sm a ll Spi tzer unstable sub- 
systems. IGiersz! (|200lD and lGierszl (|2006l ) have shown that 
the Monte Carlo code accurately re-produces binary burning 
in the cluster core, indicating that the code can resolve the 
statistical properti es of strong inter a ctions in dense regions 
well. Furthermore IHeggie fc Gierszl l|2009l ) have compared 
direct TV-Body and Monte Carlo simulations for the case of 
the cluster NGC 6397 and have shown good agreement both 
for the escape rate (driven by ejections from the core) and 
binary energy generation. Thus despite its limitations, the 
Monte Carlo code still seems to be able to produce statisti- 
cally reliable results for strongly interacting r egions. We will 
also b e able to compare our results to those of lO'Learv et al.l 
(2006), who make the explicit assumption of a Spitzer un- 
stable BH population, in order to constrain this effect. Fur- 
thermore, if there are a large number of black holes in the 
cluster any Spitzer unstable subsystem that they may form 
will be large enough that the Monte Carlo approximation 
remains valid. 

Another issue is our treatment of strong few-body in- 
teractions. O ur code curre n tly uses analyti c cross-sections 
calcul ated i nlHeggid (Il975l). iMikkolaj (1 1984 ), IStodolkiewiczl 
l|l986h . and IHeggie. Hut fc McMiUanl |l996l ) for the initiali- 



sation and outcome of binary-single and binary-binary inter- 
actions and three-body binary formation. For unequal mass 
cases these cross-sections are not certain and only allow a 
limited range of outcomes. Another problem is that lacking 
explicit orbital integration, mergers are only possible if the 
stellar radii overlap at the conclusion of an interaction, ig- 
noring the effect of close approaches during the interaction. 
Thus we will probably underestimate the number of compact 
binary mergers in our simulation. We note that ideally we 
would include a few- body integrator in t he M onte Carlo code 
as has been done bv lFregeau fc Rasiol (120071) or in the con- 
text o f a gas-Monte Carlo hybrid code by Giersz fc Spurzeml 
(2003). Such work is planned for the future but will in- 
evitably slow down the code considerably, making large pa- 
rameter studies more difficult. 
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Table 1. The initial conditions for our simulations. Column 1 
gives the model, column 2 the initial binary fraction, column 3 
the metallicity, column 4 the initial ratio of tidal to half-mass 
radius, column 5 the initial mass, and column 6 the initial half- 
mass relaxation time. Both columns 5 and 6 are averaged all ten 
independent realisations. 







Init: 


tal condit: 


ions 












Simulation 


h 


Z 


rt/r h 


M[M S ] 




[Myr] 


10sol21 


0.1 


0.02 


21 


3.61 


X 


10 5 


3.54 


X 


10 3 


10sol37 


0.1 


0.02 


37 


3.63 


X 


10 6 


1.51 


X 


10 3 


10sol75 


0.1 


0.02 


75 


3.62 


X 


10 5 


5.25 


X 


10 2 


10soll80 


0.1 


0.02 


180 


3.63 


X 


10 5 


1.41 


X 


10 2 


50sol21 


0.5 


0.02 


21 


5.08 


X 


10 5 


2.99 


X 


10 3 


50sol37 


0.5 


0.02 


37 


5.08 


X 


10 5 


1.28 


X 




50sol75 


0.5 


0.02 


75 


5.06 


X 


10 5 


4.44 


X 


10 2 


50soll80 


0.5 


0.02 


180 


5.09 


X 


10 5 


1.19 


X 


10 2 


101ow21 


0.1 


0.001 


21 


3.60 


X 


10 5 


3.55 


X 


10 3 


101ow37 


0.1 


0.001 


37 


3.62 


X 


10 5 


1.51 


X 


10 3 


101ow75 


0.1 


0.001 


75 


3.62 


X 


10 5 


5.25 


X 


10 2 


101owl80 


0.1 


0.001 


180 


3.63 


X 


10 5 


1.41 


X 


10 2 


501ow21 


0.5 


0.001 


21 


5.08 


X 


10 5 


2.99 


X 


10 3 


501ow37 


0.5 


0.001 


37 


5.07 


X 


10 5 


1.28 


X 


10 3 


501ow75 


0.5 


0.001 


75 


5.07 x 


10 5 


4.44 


X 


10 2 


501owl80 


0.5 


0.001 


180 


5.07 x 


10 5 


1.19 


X 


10 2 



3 CLUSTER MODELS 

We have preformed simulations of star clusters with 16 
different sets of initial conditions. Each simulation has 
5.0 x 10 J centres of mass(single stars or binaries). All 
simulations use a Kroupa ini tial mass function (IMF) 
l|Kroupa. Tout fe Gilmorej|l99oj ). a broken power-law with 
a low-mass slope of q ( = 1.3, a high-mass slope of a h = 
2.3 and a bre a k ma ss of Mbreak = 0.5M Q . We follow 
ISadowski et al.l (|2008T ) and choose masses between 0.1 Mq 
and 150 Mq. All simulations are initialised as Plummer 
models with a tidal cu t-off at rt c = 150 pc. According to 
the classical formula of ISpitzeil (|l987l l: 



where Mc is the mass of the cluster, Mg is the mass of the 
galaxy, and Rg is the galactocentric radius. For a galactic 
mass of « 6 x 10 10 M o and our cluster masses (Table [TJ this 
yields Rg ~ 9 — 10 kpc, a distance just beyond the solar 
orbit. Since we do not include disk shocking in our models, 
these represent halo clusters. We choose relatively isolated 
initial conditions to ensure that the effects we observe are 
due to internal cluster dynamics and yet can still investi- 
gate escapers. The tidal cut-off is not held constant dur- 
ing the evolution of the cluster but is re-calculated at each 
timestep according to the current cluster mass. There is an 
iV-dependent parameter, a, that describes a modification of 
the tidal radius necessary to produce an agreement in escape 
rate between direct iV-body models and the Monte Carlo 
code <|Baumgardtll200ll ; lGiersz. Heggie fc Hurlevll2008h . It is 
set to 1.31 in our simulatio ns which is the va l ue ch osen for 
the M4 models described in iHeggie fc Giers j (|200ct ). These 
models have a similar number of particles to the simulations 
described in this paper and agree well with observations. 

We use two different metallicities for our simulations: 
Zh = 0.02 and Zi = 0.001. Z% corresponds roughly to 
solar metallicity while Zi corresponds both to the low- 
metallicity peak of the galactic globular cluster distribu- 



Table 2. Number of BHs formed in each model by stellar evolu- 
tionary processes. N s gu is the total number of black holes formed 
in the cluster, NbBH the total number of black holes formed in 
binaries. Ngjjgu is the total number of black hole-black hole bi- 
naries formed by stellar evolutionary processes and N surv is the 
number of binaries that form a single black hole and survive the 
formation process. All primordial BH-BH binaries are disrupted 
during the second supcrnovae. Each quantity is averaged over all 
ten independent realisations and includes the rms scatter. 





Individual Black Hole Statist: 


ICS 




Simulation 


N sB H ± " 


N bBH ± CT N BHBH ± " 


N sur v ± a- 


10sol21 


1103 ± 18 


197 ± 8 


2 ± 1 


± 


10sol37 


1119 ± 44 


206 ± 17 


3 ± 2 


± 


10sol75 


1104 ± 31 


200 ±12 


2 ± 1 


± 


10soll80 


1129 ± 22 


190 ± 18 


3 ± 1 


± 


50sol21 


1495 ± 43 


976 ± 35 


13 ± 5 


± 1 


50sol37 


1515 ± 32 


989 ± 28 


12 ± 4 


± 


50sol75 


1498 ± 28 


978 ± 29 


9 ± 2 


± 1 


50soll80 


1555 ± 42 


956 ± 36 


9 ± 3 


± 


101ow21 


1248 ± 34 


217 ± 19 


2 ± 2 


± 1 


101ow37 


1262 ± 30 


228 ± 18 


4 ± 2 


± 1 


101ow75 


1265 ± 30 


229 ± 14 


5 ± 4 


± 1 


101owl80 


1296 ± 37 


227 ± 18 


3 ± 2 


± 1 


501ow21 


1719 ± 46 


1090 ± 25 


8 ± 2 


3 ± 2 


501ow37 


1728 ± 47 


1124 ± 45 


17 ±4 


3 ± 2 


501ow75 


1731 ± 37 


1125 ± 39 


16 ± 3 


2 ± 1 


501owl80 


1769 ± 31 


1077 ± 33 


9 ± 3 


4 ± 3 


With 50 Realisations 


10sol75 


1099 ± 29 


201 ± 14 


2 ± 1 


± 


501ow75 


1746 ± 38 


1132 ± 29 


18 ± 4 


3 ± 2 



tion a nd, for comparis o n pur poses, to the metallicity cho- 
sen by ISadowski et ahl ((2008). These two metallicities also 
fall within the high- and low-mass peaks in the bimodal 
metalli city distribu ti on fo und for brightest cluster galax- 
ies by lHarris et alJ l|2006h . For the purpose of our study 
the primary difference between these metallicities is the 
treatment of stellar mass-loss in the BSE stellar evolution 
code. The mass of single BHs is c alculated according to 
iBelczvnski. Kalogera fc Bulikl (|2002h . In these prescriptions 
mass-loss is suppressed at low metallicity due to less efficient 
line-driving of stellar winds and this allows the formation of 
significantly more massive BHs. These high-mass BHs mass- 
segregate more swiftly than their low-mass counterparts and 
will be stronger gravitational wave sources. The initial dis- 
tribution o f BH masses these pres criptions yield is studied 
in detail in IBelczvnski et all lj200rj ). 

We also use two different binary fractions: fb = 0.1 and 
fb = 0.5. Thus, while they have the same total number of 
centres of mass, simulations with ft = 0.1 have 5.5 x 10 5 
stars whereas simulations with ft = 0.5 have 7.5 x 10 5 
stars. The simulations with ft = 0.5 will be more mas- 
sive and produce a larger total number of BHs simply be- 
cause there are more stars present. /;, = 0.5 will also pro- 
duce more binary-single and binary-binary interactions due 
to the larger number of primordial binaries. This in turn 
will increase the probability of exchange interactions and 
increase the number of BHs introduced into binaries. The 
initial binary parameters are produc ed using t he eig envalue 
evolution and feeding algorithms of iKroupal (|l995f ). These 
prescriptions use a thermal distribution of birth eccentrici- 
ties (f(e b) = 2e&), birth mass ratios (tfc) drawn at random 
from the iKroupa. Tout fc Gilmord l|l993h IMF, and a birth 
period distribution of: 
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Time [Gyr] 

Figure 1. Run 501ow75 for two values of 7: 7 = 0.02 (solid) and 
7 = 0.11 (dotted). 7 = 0.11 has a slightly shorter relaxation time 
and produces slightly more BH-BH binaries after a Hubble time 
but both agree to within the rms error. 



f(P b ) = 2.5- 



logP(,-l 



(2) 

45 + (log Pi, -l) 2 V ' 

with the limits log P btmin = 1 and log Pb,max = 8.43. These 
birth values are then modified according to the eigenvalue 
and feeding algorithm to simulate the effect of pre-main se- 
quence evolution. Initial eccentricities are calculated as: 



lne in = — p + hie b 



where 



dt p 



(3) 



(4) 



where p'" 1 is the circularisation timescale, At « 10 5 yr is 
the pre-main-sequence evolution timescale, 7? por i is the peri- 
center distance of the binary, and A = 28 and \ — 0-75 are 
empirically determined constants. The initial mass ratio is 
given by: 



qb + (1 - Qb)p* 



where 



P<1 
p>l 



(•») 



(6) 



where the mass of the secondary is modified according to 
rri2,in ~ qi n rri2,b and the mass of the primary is unchanged 
fn\,in = mi,b. Finally the period is given by: 

/',„ = /',,[ — ) 1 f^- £ M 3/2 (7) 



/ m t ,b 

\m t ,in 



where m t j, and mt,m are the total masses before and af- 
ter the application of Equation [5] The main effect of the 
eigenvalue and feeding evolution is to depopulate the short- 



period, high-eccentricity area of the period-eccentricity di- 
agram as observed in galactic binaries. The initial period 
distribution and the eff ect of the eige nvalue feeding can be 
seen in figures 1 and 2 of lKroupal (|l995l ). The effect of a dense 
stellar evironment on generic binary p opulatio n s has been 
investigated in severa l sour ces su ch as |Heggiel 1 19751 ) and 
iPortegies Zwart et al.l (|2004l ). The lKroupal ij 19951 ) prescrip- 
tions also provide a good match to the binary parameters 
observed in the galactic field. 

We have performed simulations with four initial concen- 
trations that we control using the ratio of the initial tidal 
radius to the initial half-mass (r^). We use initial ratios of 
ft/fh — 21, 37, 75, and 180, corresponding to initial number 
densities within ru of ~ 10 2 , 10 3 , 10 4 , and 10 s respectively. 
The initial concentration prim arily affects th e half-mass re- 
laxation time (t r h), defined by Spitzcr (1987) as: 



t rh = 0.138- 



N 



1/2 3/2 



m) 1 /2GV2i n7 Ar 



(8) 



Where N is the total number of centres of mass in the sys- 
tem, (m) is the average mass, In 7 is the Coulomb loga- 
rithm, and 7 is an empirically fitted parameter. The value 
of 7 is different for single and multimass systems and has 
been debated in the literature. F or equal mass system s 
7 = 0.11 seems to be favoured l|Giersz fc Heggia Il994l ) 
whereas for multimass systems 7 = 0.02 gives better results 
ijGiersz. Heggie fc Hurlevll2008l ). We choose 7 = 0.02 for our 
simulations but have re-run one set, 501ow75, with 7 = 0.11 
for comparison purposes. The results will be discussed in §|4] 
The simulations all have a fixed initial tidal cut-off and thus 
the more concentrated simulation have smaller values of r^. 
Thus, according to Equation |8j they will also have shorter 
values of t T h and will evolve, dynamically speaking, more 
quickly than their less concentrated counterparts. Whatever 
effect dynamics have on the production of BH-BH binaries 
should be accelerated in these systems. 

Since star cluster dynamics are stochastic and chaotic, 
large fluctuations can occur between different realisations 
of the same simulation (e.g. Giersz. Heggie fc Hurley! [20081 : 
iHeggie fc Gierszll2008l : iGiersz fc Heggiell2009D . For this rea- 
son we perform ten independent realisations of each com- 
bination of initial conditions differing only by the initial 
random seed. Thus we have a total of 160 simulations to 
analyse. To ensure that ten simulations is enough for con- 
vergence we have extended two of the simulations, 10sol75 
and 501ow75, to 50 simulations and will compare the num- 
ber of BH-BH binaries produced in § [4] Table [T] gives the 
initial parameters of our 16 different sets of simulations, av- 
eraged over all ten realisations. Each simulation is run on a 
single processor at the HLRS supercomputer in Stuttgart. 
The shortest simulations (10sol21) take ~ 4 h to complete 
and the longest (501owl80) take ~ 12 - 16 h. 



4 RESULTS 

We find no more than one or two NS-NS or NS-BH binaries 
over the course of an entire Hubble time in any of our simu- 
lations and no NS-NS or NS-BH mergers. For this reason we 
only present results for the BH-BH binary population. The 
lack of NS-NS and NS-BH binaries is due to the fact that few 
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Table 3. The cumulative number of BH-BH binaries after 3, 9, 
and 25 t r h, and also after ITh- Each column is averaged over all 
ten independent realisations and includes the rms scatter. A dash 
in a column indicates that the cluster did not reach that number 
of half-mass relaxation times within one Hubble time. 





BH-BH ] 


3inarics Aft> 






Simulation 


t = 3t rh 


t = 9t rh 


t = 25t rh 


t = 14 Gyr 


10sol21 


1 ± 1 






1 ± 1 


10sol37 


1 ± 1 


14 ± 11 




14 ± 11 


10aol75 


± 1 


8 ± 6 


49 ± 19 


52 ± 19 


10soll80 


± 


12 ± 6 


54 ± 21 


123 ± 27 


50sol21 


1 ± 1 






3 ± 2 


50sol37 


3 ± 2 


36 ± 10 




50 ± 11 


50sol75 


1 ± 1 


26 ± 8 


115 ± 24 


147 ± 28 


50soll80 


± 1 


11 ± 5 


111 ± 23 


354 ± 33 


101ow21 


22 ± 10 






27 ± 10 


101ow37 


23 ± 4 


44 ± 6 




44 ± 6 


101ow75 


18 ± 10 


32 ± 13 


54 ± 6 


54 ± 16 


10Iowl80 


26 ± 8 


51 ± 8 


79 ± 20 


112 ± 24 


501ow21 


104 ± 16 






127 ± 16 


501ow37 


93 ± 22 


175 ± 29 




184 ± 29 


501ow75 


64 ± 10 


155 ± 22 


173 ± 22 


202 ± 38 


501owl80 


103 ± 19 


205 ± 35 


294 ± 50 


453 ± 109 




With 


50 Realisat 


ions 




10aol75 


± 1 


7 ± 4 


38 ± 20 


41 ± 21 


501ow75 


78 ± 15 


175 ± 29 


232 ± 43 


245 ± 48 



primordial binaries survive to a phase where they would con- 
tain an NS or BH since most merge during mass transfer and 
those that do are disrupted at the second supernova. As we 
will show, our BH-BH binaries are not primordial but rather 
form dynamically, a process that occurs most efficiently for 
massive objects. Since the BHs are, for the most part, signif- 
icantly more massive than the NSs in our simulations, they 
will be preferentially exchanged into binaries until they are 
depleted. None of our simulations are completely depleted 
of BHs after one Hubble time and thus the NSs have little 
opportunity to take part in dynamical compact binary for- 
mation. We do not analyse the (large) WD- WD population 
but save these results for a future paper. 

In Table [2] we show the total number of BHs formed 
both in isolation and in binaries in each simulation aver- 
aged over all ten independent realisations. For each simu- 
lation the rms scatter across the independent realisations 
is small and is merely a result of random sampling of the 
IMF. The number of BHs formed is not a function of con- 
centration in any range of Z or /;,. This is not surprising 
since BHs are produced primarily by individual stellar evo- 
lutionary processes. It is possible that extra BHs could be 
produced by stellar collisions in which two stars below the 
critical mass to produce a BH merge to form a single star 
above the critical mass. Collisions would be expected to be 
more frequent in dense clusters but this effect does not seem 
to produce a significant enhancement in the number of BHs 
in our simulations. The number of BHs formed depends on 
fb because a higher binary fraction corresponds to a larger 
number of stars, and on Z because mass-loss is less efficient 
at low metallicity and stars with a lower zero age main se- 
quence mass can become BHs. Proportionately more BHs 
are formed in binaries at = 0.5 than at fb = 0.1 but this 
is simply a consequence of the larger fraction of stars found 
in binaries at high fb. It is apparent from column 5 of ta- 
ble [2] that very few of the binaries that form a single BH 



survive is formation; most either merge or are disrupted at 
the supernovae. Furthermore, very few binaries where both 
members are black holes form directly from primordial bi- 
naries and of those all are disrupted during the formation of 
the second BH. Therefore all BH-BH binaries produced by 
our simulations must be formed by dynamical means. There 
is no difference in either the mean number or the rms scat- 
ter in BHs produced in either of the simulations for which 
we have performed an extended number of realisations. We 
therefore conclude that ten simulations are sufficient to pro- 
duce accurate statistics on the BH population produced by 
stellar evolution. 

In Table [3] we present the cumulative number of BH- 
BH binaries that have existed in each simulation up to the 
dynamical time shown and after one Hubble time (Th = 14 
Gyr). The ages 3t r h, 9t r h, and 25t r h correspond respectively 
to the dynamical age of the clusters with rt/r^ = 21, 37, and 
75 at ~ 1Tb- The simulations with rt/ru = 180 have a dy- 
namical age of ~ 115— 120t r h at ITh- A new BH-BH binary 
is counted every time a binary where both members are BHs 
forms. Both a binary with one BH and one main sequence 
(MS) star where the MS star is exchanged for a BH and a 
binary where both members are BHs and one of the BHs is 
exchanged for a new BH are counted as new BH-BH binaries. 
The rms scatter in the number of BH-BH binaries produced 
by different realisations of the same simulation in Table [3] is 
significantly larger than for Table [5] This is because, unlike 
stellar evolution, dynamical binary formation is a stochas- 
tic process, strongly dependent on chance encounters that 
are a function of the detailed dynamics of the specific sys- 
tem, and, since all BH-BH binaries are formed dynamically, 
large system-to-system variations are expected. Again we 
find that for the two simulations where we performed ad- 
ditional independent realisations, the additional realisations 
provide no difference in the size of the rms scatter and the 
average number of BH-BH binaries at each relaxation time 
is the same to within the rms scatter. Therefore we conclude 
that ten realisations is also enough to constrain the statis- 
tics of the dynamical processes and we retain this number 
for the rest of our analysis. 

There are some clear trends in Table [3] After ITh the 
number of BH-BH binaries increases with and initial con- 
centration and decreases with Z. The reasons for the cor- 
relation are clear: a larger number of both BHs and hard 
binaries for the BHs to be exchanged into. The correlation 
with initial concentration is related to the relative dynami- 
cal ages of the clusters after ITh- For the same fb and Z the 
simulations have roughly the same number of BH-BHs when 
compared at the same dynamical age. By ITh, however, 
clusters with a higher initial concentration are dynamically 
older and have had more opportunity to produce BH-BH 
binaries than their dynamically younger counterparts. The 
correlation with Z is due to the higher mass BHs present 
at low metallicity. Since the mass-segreg ation timescal e, t aq , 
(for a two-component system) scales as (|Spitzer|| 19871 ): 

teq OC t rh — (9) 

m2 

where rrt2 > mi llWatters. Joshi fc Rasid l2000l : 
iKhalisi. Amaro-Seoane fc Spurzem! l2007h . we expect 
more massive BHs to mass segregate more rapidly than less 
massive ones. This accelerates the process of BH-BH binary 
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Figure 2. The number of BH-BH binaries per Gyr for each of the 16 sets of initial conditions. Each time bin is averaged over all ten 
independent realisations and the error bars represent the rms scatter. 



formation since the core is the densest region and is where 
most dynamical binary formation will take place. The more 
massive BHs are also more likely to be retained by the 
cluster since they will need larger kicks upon formation in 
order to reach escape velocity. 

We have also considered a different value of 7 and we 
present these results in figure[TJ According to equation^ the 
effect of increasing 7 is to reduce t T h by the logarithm of the 
same factor. Thus the simulations with 7 = 0.11 are slightly 
more dynamically evolved after ITh that the simulations 
with 7 = 0.02. Consequently these simulations produce a 



few more BH-BH binaries. The increase, however, is small 
and the average values for both sets of simulations fall within 
each others rms scatter. We conclude that any reasonable 
changes in 7 will not affect our results in any important 
way. 

Figure [2] shows the number of BH-BH binaries in each 
simulation per Gyr. The trends noted in Table [3] are appar- 
ent, particularly those associated with the initial concentra- 
tion. Simulations with higher initial concentration have a 
peak number of BH-BH binaries per unit time much ear- 
lier than those with lower concentration. Those with lower 
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Figure 3. The half-mass radii for the BH system. Top to bottom: Z = 0.02 and f b = 0.1, Z = 0.02 and f b = 0.5, Z = 0.001 and f b = 0.1, 
and Z = 0.001 and f b = 0.5. Left to right: rt/r^ = 21, 37, 75 and 180. Shown are cluster (dotted), of all BHs (dashed), and 
of all BH-BH binaries (solid). The radial profiles have been boxcar smoothed with a box 100 Myr box and each box has been averaged 
across all ten independent realisations. 



concentration sustain more constant but lower BH-BH pop- 
ulations over longer spans of physical time. It is also clear 
that the more metal-poor simulations evolve more quickly 
and produce BH-BH binaries earlier than their metal-rich 
counterparts. This again is a consequence of faster mass seg- 
regation. Finally, ft does not affect the time of peak BH-BH 
binary number but simulations with = 0.5 are able to 
sustain production of BH-BH binaries longer because the 
supply of BHs and hard binaries to exchange them into is 
larger in these clusters. 



In Figure[3]we investigate the spatial distribution of the 
BHs and BH-BH binaries in each simulation by comparing 
the half-mass radius of each species to the half-mass radius 
of the entire cluster. In all cases single BHs are centrally con- 
centrated compared to the rest of the stars. This is simply 
a consequence of mass segregation. The varying mass segre- 
gation timescales can be seen in the length of time taken for 
the half-mass radii of the BH populations to contract to an 
equilibrium state caused by binary burning. The half-mass 
radii of the BH-BH populations must be interpreted more 
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Figure 4. The cumulative chirp mass distribution of BH-BH binaries up to ITu binned in IMq bins. Each bin in each panel is averaged 
over all ten independent realisations and the error bars give the rms scatter. 



carefully. Since BH-BH binaries are formed dynamically in 
cluster cores, the BH-BH binary populations is initially very 
centrally concentrated. As the population evolves, however, 
the BH-BH binaries interact strongly and can eject each 
other from the core region. As is clear from Figure [2] there 
are often very few BH-BH binaries in the cluster, even over 
a span of 1 Gyr, and thus the half-mass radius of the BH- 
BH population in Figure [3] is often based on quite a small 
number of objects. Thus the location of a single massive BH- 
BH binary can dominate the determination of the half-mass 
radius of the BH-BH binary population. Overall it is clear 



that BH-BH binaries form in the cluster core and tend to 
be centrally concentrated but they are not necessarily more 
concentrated than the single BH population. Individual BH- 
BH binaries may also exit in the outskirts of the cluster for 
a time if they are ejected from the core by dynamical inter- 
actions and before they have a chance to sink back to the 
center due to dynamical friction. 

It is also worth noting that the BH-BH binaries will 
not be strongly affected by interactions with anything other 
than BHs. The comparative masses mean that other stars do 
not have sufficient energy to disrupt or scatter the BH-BH 
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Figure 5. The cumulative chirp mass distribution of BH-BH binaries up to 9t r h for six simulations binned in IMq bins. Both plots are 
for Z = 0.001 with the left plot having fi, = 0.1 and the right plot having = 0.5. Concentrations are rt/rh = 37 (solid), rt/rh = 75 
(dotted), and rt/rh = 180 (dashed). 9t r h is the dynamical age of the simulation with rt/rh = 37 after one Hubble time. The simulations 
with rt/rh = 21 are not shown since they do not reach 9*, ^ within one Hubble time and have too few BH-BHs at 3t r ^ for interesting 
statistics. Each bin is averaged over all ten independent realisations. 



binaries. Conversely the BH-BH binaries will have a very 
strong effect on the other stars they encounter and will be 
a the major energy source in the cluster core. Thus BH-BH 
binaries are only affected by each other but have a major 
influence on core dynamics. 

The cumulative chirp mass (M c hi rp = 
(M 1 M 2 ) 3/5 /(M 1 + M 2 ) 1/5 ) distribution for all BH-BH 
binaries that have existed in the simulations up to ITh 
are given in Figure [4] We choose to display M cri i rp rather 
than the total mass because the amplitude of gravitational 
radiation, ho, scales as 



ho oc 



G 5/3 w2 /3 M S/3 



chirp 



(10) 



|Pierro et al.ll20o"lh . Thus it is M c hi rp and not the total mass 
of the binary that is significant for gravitational wave detec- 
tion. For the high metallicity simulations the distribution of 
chirp masses is narrow with a peak around 8 — IOMq . The 
distribution is not affected by any of the other initial condi- 
tions. The low-metallicity distribution is broader and fairly 
flat between 10 — 20Mq . This is a direct consequence of the 
more massive BHs generated at lower metallicity. Here the 
distribution is weakly affected by the initial concentration 
with M c hirp peaking at lower masses for the more concen- 
trated simulations. This is a result of the relative dynami- 
cal ages of the simulations as we demonstrate in Figure O 
the distribution of M crl ir P after 9t r h- At the same dynamical 
age the mass distribution is unaffected by the concentra- 
tion. Recalling Equation [9] the more massive BHs will mass 
segregate before the less massive BHs and thus interact and 



be disrupted or ejected earlier. Only then will low-mass BHs 
participate in BH-BH binary formation. Since the more con- 
centrated clusters are dynamically older, they have had more 
time to experience this effect, deplete their high-mass BHs, 
and have a lower-mass BH-BH binary population. The high- 
metallicity clusters do not have a broad enough distribution 
in mass for this effect to be important. Perhaps the most 
interesting result is that after ITh the M c hi rp distributions 
are systematically different between clusters with different 
metallicities and concentrations. Building an M cri i rp distri- 
bution from gravitational wave observations can yield infor- 
mation on the physical and dynamical state of GCs. 

The cumulative distribution of BH-BH binary binding 
energy (-Ebind) up to ITh is given in Figure [6] The energy 
is given in units of the thermal energy of the core of the 
cluster, 



k B T = (m CO rc/2A r core)cr c 



(11) 



where m corc is the total core mass, N cole the number of stars 
in the core and <r C orc the velocity dispersion in the core. The 
core quantities are chosen because the BH-BH binaries are 
formed and interact in the core, making this region the most 
relevant for the dynamics. All BH-BH binaries are hard. 
This is to be expected because soft binaries would be de- 
stroyed by the interactions necessary to introduce a BH into 
them. There is little variation with the cluster parameters 
because ksT, the normalisation factor for -Ebind, scales with 
core mass and density. The only exception to this appears 
to be clusters with rt/rh = 180 where there are an excess 
of soft binaries. The combination of the larger interaction 
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Figure 6. Cumulative binding energy distribution of BH-BH binaries up to lTu binned uniformly in log space. Each bin is averaged 
over all ten independent realisations and the error bars give the rms scatter. Energies are in units of the thermal energy of the core of 
the cluster. 



cross-section for binaries with larger semi-major axes and 
the larger interaction rate in the cores of very dense clusters 
may lead to more BHs being exchanged into softer bina- 
ries in the high-density simulations. It could also be that 
the soft binaries in the more concentrated simulations still 
have a larger binding energy in physical units than those in 
the less concentrated simulations and have a slightly better 
chance of survival. 

In Figure [7] we present the cumulative period (P) dis- 
tribution of all BH-BH binaries up to 1Tb ■ There is a large 
spread in P, ranging from days to ~ 10 5 yrs. The period 



distribution does not depend strongly on the cluster pa- 
rameters. The entire distribution is shifted slightly towards 
short periods in the more concentrated simulations. This 
is partly due to the slightly higher velocity dispersion in 
these clusters and the consequently higher value of fcgT in 
physical units. Thus in concentrated clusters binaries must 
have higher binding energies and consequently shorter peri- 
ods in order to be above the hard-soft boundary (such as it 
is in multi-mass systems). It is also partly a product of the 
larger number of hardening interactions due to the higher 
interaction rate. Despite the overall shift, the clusters with 



12 J. M. B. Downing, M. J. Benacquista, M. Giersz, & R. Spurzem 



-2 024 -2 024 -2 024 -2 024 



2 0.1 

z 

\ 

0.01 



2 0.1 

z 

0.01 



2 0.1 

\ 

go 0.01 



3 0.1 r 



2 

^ 0.01 




-2 2 4 
log(P) [Yrs] 



-2 024 -2 024 -2 024 
log(P) [Yrs] log(P) [Yrs] log(P) [Yrs] 



Figure 7. Cumulative period distribution of BH-BH binaries up to ITjj binned uniformly in log space. Each bin is averaged over all ten 
independent realisations and the error bars give the rms scatter. 



Tt/rn = 180 show a peak in at the longer end of the distri- 
bution. This simply reflects the peak at low binding energy 
seen in Figure [6] Overall, however, the period distributions 
are similar and span approximately the same range for all 
models. 

Although there are some binaries with periods less than 
a year present in most simulations, most binaries do not 
have sufficiently short periods to p roduce g ravita tional wave 
mergers within ITh- According to iPetersl (|l964l ) the rate of 
semi-major axis decay, a, for a binary emitting gravitational 
waves in the orbit-averaged regime is: 



<d) = 



64 G raim.2 (mi + 012) 
5 C 5a3(f-e 2 ) 7/2 



. 73 2 37 4 > 

1 H e H e 

74 96 1 



(12) 



where mi is the mass of the primary, 7712 is the mass of 
the secondary, e is the eccentricity, and a is the semi-major 
axis. We can calculate an approximate inspiral time, tin, f° r 
a binary by taking: 



fin 



(13) 



where Ctin IS the initial semi-major axis of a binary. For a 
circular binary with mi = 7712 = WMq and an initial period 
of P in = 1 day, t in w 1 Gyr. If P ln = 1 yr then t- m = 10 6 
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Gyr. Thus all but the shortest period binaries in our clus- 
ters will be unable to merge within ITh- Furthermore, the 
minimum gravitational wave frequency for which LISA is 
sensitive is « 10 -5 Hz. For circular bi naries all power is 
emitted in the n = 2 harmonic of u> (IPeters fc Mathews! 
ll963l ; lBelczvnski, Benacquista fc Bulildl2010h . Thus for cir- 
cular binaries to be detected by LISA they must have periods 
of less than a day. 

The presence of eccentricity in a binary can vastly 
improve its prospect fo r gra vitational wa ve detection. 
IPeters fc Mathews! l|l963h and iPeters! (QUI) have shown 
that with increasing eccentricity gravitational waves are 
preferentially emitted at higher harmonics of the orbital 



frequency. This can enhance the power emitted in gravi- 
tational waves by a factor of 10 2 or more for binaries with 
e > 0.8 and reduce ti n by a similar factor. This enhances 
the chance of a relativistic merger within ITh- Emission at 
higher orbital harmonics produces gravitational waves with 
higher frequencies that would be produced by identical bi- 
naries with circular orbits. Thus eccentricity can bring long- 
period binaries into the LISA band. In Figure [8] we show 
the eccentricity of our BH-BH binaries as a function of pe- 
riod. It must be noted that these eccentricities are produced 
randomly upon creation of a binary, do not experience self- 
consistent dynamical evolution during interactions, and are 
subject to circularisation in the course of binary evolution in 
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Figure 9. Cumulative distribution of BH-BH binary gravitational wave inspiral timescale up to ITh binned uniformly in 
Each bin is averaged over all ten independent realisations and the error bars give the rms scatter. 



space. 



BSE. Interactions tend to increase the eccentricity of bina- 
ries and thus the eccentricities produced by our simulations 
should be taken as lower limits. Even with this caveat, Fig- 
ure [8] shows that there are a wide range of eccentricities for 
any given period. 

Since some of our binaries are eccentric and since this 
eccentricity can significantly reduce the inspiral timescale 
of the binary, we use Equation [13] to estimate the inspiral 
timescale of all BH-BH binaries in our simulations. The re- 
sult is given in Figure [9] the trends of which simply reflect 
the period distribution in Figure [7] It is apparent that even 



with eccentricity included there are very few binaries with an 
inspiral timescale shorter than ITh- For the few BH-BH bi- 
naries with < 10 5 Myr the dynamics play a destructive 
role. Figure [10] shows the timescale for dynamical disrup- 
tion or ejection of BH-BH binaries, t<ji S . It is apparent that 
the average tdis for BH-BH binaries is very short, between 
1 — 100 Myr in most cases, and is shorter in the more con- 
centrated clusters due to the higher interaction rate. There 
is little between the distributions in Figure [5] and Figure [TO] 
and ti n > tdis in almost all cases. Indeed in none of our 
simulations do we find any BH-BH binary mergers. Some of 
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Figure 10. Cumulative distribution of BH-BH binary dynamical disruption timescale up to ITh binned uniformly in log space. Each 
bin is averaged over all ten independent realisations and the error bars give the rms scatter. 



the eccentric, short-period binaries in our simulations may, 
however, have a chance of appearing in the LISA band and 
we will turn to this possibility in § 15.21 

Although we find no mergers within the clusters, the 
simulations eject hard binaries. Figure [TTI shows the distribu- 
tion of binding energies for escaping BH-BH binaries. i?bmd 
is, on average, much higher for the escapers than for the 
system as a whole. This is because the most tightly bound 
binaries tend to receive the highest recoil velocities in few- 
body encounters and are thus the most likely to be ejected 
from the system. Therefore the most promising merger can- 



didates are the least likely to remain in the cluster. The dis- 
tribu tion in Figure II II com pares favourably with the results 
from lO'Learv et alJ (|2006l ) (their Figure 6) which were de- 
rived based on explicit few-body integration. This agreement 
increases our confidence that both the few-body encounters 
and the dense core dynamics are being re-produced success- 
fully in our code. 
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Figure 11. The binding energy of BH-BH binary escapers binned 
uniformly in log space. The top row shows the Z = 0.02 simula- 
tions and the bottom the Z = 0.001 ones. In the right column are 
simulations with /{, = 0.1 and in the left simulations with /(, = 0.5 
column. Concentrations are rt/rf, = 21 (solid), rt/r^ = 37 (dot- 
ted), rt/r'h = 75 (dashed) and rt/r^ = 180 (dot-dashed). Each 
bin is averaged over all ten independent realisations. 



5 GRAVITATIONAL WAVE DETECTION 

In this section we consider the prospects for gravitational 
wave detection in both the ground- and space-based fre- 
quency bands. 



5.1 Ground-Based Sources 

Due to the long inspiral times and short disruption times 
shown in figures and [TO] we find no compact mergers in 
any of our simulations. Thus the BH-BH binary destruction 
rate within the cluster dominates the creation rate and all 
BH-BH binaries are either ejected or disrupted before they 
have a chance to inspiral and merge. We therefore predict 
that globular clusters as objects will not contribute signifi- 
cantly to the ground-based gravitational wave detector sig- 
nal. Many of the escapers produced by GC dynamics should, 
however, merge in the galactic field within a Hubble time 
and thus may contribute to the detection rate. We save anal- 
ysis of these sources for a follow-up paper and concentrate 
instead on the very interesting LISA sources that we find 
within our simulations. 



5.2 LISA Sources 

In order to study potential LISA sources from these simula- 
tions, we consider only those binaries present in the cluster 
simulations at ages between 10 and 14 Gyr. This is covers 
the age distribution of GCs in the Milky Way and stellar 
mass inspirals will not be strong enough gravitational wave 



sources to be observed beyond the galaxy. Advanced LIGO 
will be able to detect BH-BH inspirals out to moderate red- 
shifts and for these predictions we will have to take the ear- 
lier stages of cluster evolution and the possibility of younger 
GCs in other galaxies into account. Of these Milky Way bi- 
naries, we further restrict them to have a combination of 
orbital periods and eccentricities such that the harmonic 
with peak power lies within the LISA sensitivity band. We 
estimate the frequency of this harmonic by approximating 
the orbit at periastron (o(l — e)) with a circular orbit of 
radius r = a(l — e). The frequency of this circular orbit is 
then proportional to (1 — e) -1,5 . However, since the orbital 
speed of the eccentric binary at periapsis is higher than that 
of a circular binary, the dominant frequency of the gravita- 
tional wave burst will be somewhat higher than the circular- 
frequency. For a parabolic encounter, this speed difference 
is y/2, and so the harmonic with peak power is estimated to 
be 



(i-er 



(14) 



With this age and frequency/eccentricity restriction, we find 
that 33 simulations have potential LISA sources. All but one 
of these have only one potential LISA source at any time 
during the 10 to 14 Gyr under consideration. The one excep- 
tion (a realization of the 501ow37 model) has two potential 
LISA sources during the age span of 10.6 to 11.6 Gyr. Six 
of these simulations had binaries with e < 0.9, with one of 
these having e < 0.7. Although there are 34 potential LISA 
sources in these simulations, we need to determine if these 
BH-BH binaries will have sufficient strength to be detected 
by LISA. This is done by calculating the signal-to-noise ratio 
using: 



p 2 = 4 



\kf)i 

Sn(f) 



df 



(15) 



where h(f) is the Fourier transform of the response of LISA 
to the gravitational wave and S n (f) is the power spectral 
density of the expected noise in LISA. We include both in- 
strument noise and an estimate of the Gala ctic white dwarf 
binary foreground from iRuiter et all (|2007l ). 

We determine the response, h(t) by placing ten reali- 
sations of each binary in Galactic globular clusters that are 
within 5 kpc of t he Earth. The properties of these clusters 
are obtained from lHarrisI (Il996l ) , and are shown in Table [4] 
Each realisation is given a sky location within the globular 
cluster and then assigned random orientations and initial 
orbital phases. The barycenter ed waveform is det ermined 
using the harmonic expansion of lPierro" et al.l (|200ll ) carried 
out to the n ~ 1300 harmonic. The response of LISA is 
calculated using t he lon g wavelength approximation as de- 
scribed in lCutlerl (|199a ) for one year of observation. 

Setting a detection threshold of p 7 in a single inter- 
ferometer (which corresponds to a combined p 10 in two 
channels of the LISA data stream) we find potentially de- 
tectable binaries in three realisations of simulation 501ow37 
and one realisation of simulation 50sol75. The properties of 
these binaries are shown in Table [5] The notable features of 
these binaries are that they are highly eccentric. The long- 
period binary from realisation nine of simulation 501ow37 is 
only visible from a few orientations in the nearby globular 



Table 4. Celestial coordinates and distances for the 16 globular 
clusters within 5 kpc of Earth 



dist 
[kpc] 



NGC104 


00 


24 


05.2 


-72 


04 


51 


4.5 


E3 


09 


20 


59.3 


-77 


16 


57 


4.3 


NGC3201 


10 


17 


36.8 


-46 


24 


40 


5.0 


NGC6121 


16 


23 


35.5 


-26 


31 


31 


2.2 


NGC6218 


16 


47 


14.5 


-01 


56 


52 


4.11 


NGC6254 


16 


57 


08.9 


-04 


05 


58 


4.4 


NGC6366 


17 


27 


44.3 


-05 


04 


36 


3.6 


NGC6397 


17 


40 


41.3 


-53 


40 


25 


2.3 


NGC6540 


IS 


06 


08.6 


-27 


45 


55 


3.7 


NGC6544 


18 


07 


20.6 


-24 


59 


51 


2.7 


2MSGC01 


18 


08 


21.8 


-19 


49 


47 


3.6 


2MSGC02 


18 


09 


36.5 


-20 


40 


44 


4.0 


Tcrl2 


18 


12 


15.8 


-22 


44 


31 


4.8 


NGC6656 


18 


36 


24.2 


-23 


54 


12 


3.2 


NGC6752 


19 


10 


52.0 


-59 


59 


05 


4.0 


NGC6838 


19 


53 


46.1 


18 


46 


42 


4.0 



Table 5. Properties of the detectable LISA sources from these 
simulations. If the orbital period decreases during the observation 
time, a range is given. 



Simulation 


Age 


Ml 


M 2 




e 




[Gyr] 


[M ] 


[M ] 


[X10 3 s] 




501ow37.4 


< 10.2 


20.5 


25.9 


920-490 


0.988 


501ow37.9 


> 12.5 


11.1 


16.8 


1190 


0.986 


501ow37.10 


> 11.6 


28.4 


14.3 


215 


0.933 


50sol75 


< 10.1 


10.8 


8.5 


479-257 


0.998 



cluster NGC6121, while the binaries in the other realisations 
of simulations 501ow37 and 50sol75 are visible in all of the 
globular clusters chosen. The spectra of the binaries in reali- 
sations nine and four for simulation 501ow37 from NGC 6121 
are shown in figure [12] along with an estimate of the com- 
bined instrument and Galactic white dwarf binary confusion 
noise. It is also interesting to note that there is a preference 
for high primordial binary fractions, but otherwise there is 
no clear dependence on cluster parameters for these sources. 
It must again be noted that these eccentricities are not self- 
consistently produced by few-body calculations and must 
be interpreted with caution. Furthermore, the simulations 
indicate at most two binaries per globular cluster with the 
potential for detection by LISA and the bulk of the simu- 
lations result in no detectable systems. Consequently it is 
difficult to make any firm predictions about the likelihood 
of detection of BH-BH systems in the Galactic globular clus- 
ter system with LISA. Nonetheless, if we assume that high 
binary-fractions are common in globular clusters and that 
the initial concentrations are rt/ru = 37, then roughly 30% 
of nearby globular clusters may house a detectable binary. 



6 DISCUSSION 

Although we produce no compact binary mergers within 
our simulations, we can still compare o u r dyn amics to 



lO'Learv et all (120061) and ISadowski et all (|2008l ). As de- 
scribed in § [TJ lO'Learv et al.l (|2006l ) assume that the BHs 




are sufficiently mass-segregated that they form a decou- 
pled subsystem at the centre of the cluster and inte ract 
only with themselves. By contrast ISadowski et al.l (|2008l ) as- 



f(Hz) 

Figure 12. Gravitational wave spectra of two BHBH binaries 
from simulation 501ow37 compared with combined instrumental 
and Galactic white dwarf binary confusion noise. The weak sig- 
nal is from realisation 9 and the strong signal is from realisation 
4. Here, we have calculated the waveform out to the n ~ 5000 
harmonic to show the full spectrum of the binaries. For the cal- 
culation of the signal-to-noise ratio, both signals are cut-off at a 
frequencies around 2 mHz because we only use harmonics below 
n = 1300. Note that the true signal-to-noise ratio is likely higher 
than we have calculated. 



sume that the BHs always remain in equilibrium with the 
rest of the cluster and, since they do not become centrally 
concentrated, remain at much lower density than in the 
mass-segregated case . Therefore the BH-BH binaries in the 
lO'Learv et al.l (|2006h simulations frequentl y interact with 
each o ther whereas the BH-BH binaries in Sadowski ct al. 
(|2008h do not. T his means that the BH -BH binary formation 
rate reported in lO'Learv et al . (2006) is much higher than 
that reported in Sadowski et al.l ( 20081 ) . However, since the 
BH-BH binaries are the only objects in the system mas- 
sive enough to disrupt and eject other BH- BH binaries, the 
BH-BH binaries in the lO'Learv et all (2006) simulations are 
much more likely to be disrupted than the binaries in the 
I Sadowski et al.l (|2008l ) s imulations. In pr a ctice the lower dis- 
ruption rate wins and I Sadowski et alj l|2008l ) produces a 
larger and more cons tant rate of BH-BH mergers than does 
lO'Learv etal] (|200rj ). Therefore lb 'Learv et all (|2006l ) pro- 
vides a lower limit on gravitational wave detection rates and 
ISadowski et al.l (|2008l ) provides an upper limit. 

Our sim ulations, includ i ng fu ll global dynamics, sug- 
gest that the lO'Learv et al.l (|2006l ) approximation is more 
accurate. Figure [3] shows that the BHs in our simulations 
mass-segregate swiftly and that BH-BH binaries form in 
the centre of the cluster. The short dynamical disruption 
timescales in Figure [10] imply that the BH-BH binaries in- 
teract strongly with each other and are swiftly destroyed. 
The binding energy distribution of escapers in Figure QT] 
shows that hard binaries are preferentially ejected from the 
cluster, suggesting a dynamical ejection scenario. Finally 
the time-dependant number of BH-BH binaries per Gyr re- 
ported in Figure [2] particularly for the dense clusters, is 
more consistent with the time-dependent merger rates of 
lO'Learv etahl (120061) than the constant merger rate found 



Sadowski et ahl (2008) 



A major difference between our results and those of 
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both lO'Learv et al l (|2006l ) and lSadowski et alj (|2008T ) is that 
we find no BH-BH mergers within our simulated clusters. 
lO'Learv et al.l (|2006h fi nd at least 30% of thei r mergers occur 



within the clusters and Sadowski et al. (2008 ) find 90%. The 



differe nce between lo'Learv et alj 1 20061 ) and 



Sadowski et al.l 



(2008) is due to the larger number of i nteractions and hence 
the larger number of ejections in the lO'Learv et ail (|2006l ) 
simulations. Our lack of mergers within the clusters is a 
result our high interaction rates but is also affected by the 
more approximate treatment of few-body interactions in our 
Monte Carlo code. In particular our binary-bin a ry int erac- 
tion prescriptions follow those of [Stodofkicwicz (1986]) and 
in these prescriptions 88% of all interactions result in the 
disruption of the softer binary. This works well for normal 
binary-binary interactions but may over-predict the disrup- 
tion rate in interactions between two hard BH-BH binaries. 
In reality more of these binaries should probably survive and 
be hardened by the interaction. The prescriptions also gives 
0.516(i5(,i + -Em), where Ebi,2 are the binding energies of 
the two binaries, as increased binding energy to the hard 
binary and distributes the same amount as kinetic energy 
between the centres of mass after the interaction. In general 
this produces the correct hardening of the surviving binary 
(compare the energy distribution in Figure [TT] to Figure 6 in 
lO'Learv et al.l (2006)) but could well produce recoil veloci- 
ties and hence escape rates that are systematically to high. 
Furthermore these prescriptions allow neither mergers dur- 
ing the interaction nor long-lived hierarchical triples. The 
combination of all these effects means that the number of 
mergers within our clusters is almost certainly too low. 

Our simulations represent the first quantitative study 
of stellar-mass BH-BH binaries as LISA sources within star 
clusters. We can compare our results to the number of 
stellar- mass sources predicted in the galactic field popula- 
tion bv lBelczvnski. Benacquista fc Bulikl (|2010h and quoted 
in § [T] When normalised to the total mass in their simula- 
tions, iBelczvnski. Benacquista fc Bulikl (|2010h find, in the 
case where mergers on the Hertzsprung gap are not al- 
lowed, 2 x 10 -5 resolvable LISA sources in the Milky Way 
per 10 5 M Q , 1 x 10~ 5 of which per 10 5 M Q are BH-BH. In 
the case where mergers during the Hertzsprung gap are al- 
lowed they find 4 x 10" 6 per 10 5 M Q , of which none are 
BH-BH. Assuming all four of the sources in our simulations 
are resolved at the current time this yields 6 x 10 -3 de- 
tections per 1O 5 M0 from our simulations, all BH-BH. This 
rate depends on the parameters of the individual binaries 
but has no clear dependence on the cluster parameters. All 
sources are highly eccentric, indicating that eccentric stel- 
lar mass sources may be present in the LISA band. Many 
of our escaping binaries also have periods of a day or less 
and will almost certainly appear in the LISA band at some 
point. This does, however , place a limitation on the inter- 
pretat ion of the results of IBelczvnski. Benacquista fc Bulikl 
(|201Gf ). In their model the presence of stellar- mass BH-BH 
binaries in the LISA band would indicate that mergers in 
the Hertzsprung gap are rare whereas our results show that 
even if mergers on the Hertzsprung gap are common, BH-BH 
binaries could still exist in the LISA band due to dynam- 
ical processes in star clusters. However, the BH-BH bina- 
ries generated through dynamical processes will be associ- 
ated with individual globular clusters if they are retained, 
or in the halo if they have been ejected. This is a distinct 



population from the disk population o f binaries found in 
IBelczvnski. Benacquista fe Bulikl |2010l ). 

Finally we note in passing that although we have anal- 
ysed our simulations in terms of compact binaries they are 
in no way limited to such studies. Stellar evolution is cal- 
culated for all stars in the cluster and thus our simulations 
represent a complete star cluster popul ations synthesis study 
based on th e stellar evolution tracks oflHurlev. Pols fc Tout! 
(|2000l ) and iHurlev. Tout fc Pols! (|2002h . We hope to make 
the full results of our simulations publicly available in the 
near future and encourage anyone who is interested in such 
data to contact the authors. 



7 CONCLUSIONS 

We have studied the dynamics of the BH population in 
star clusters with a self-consistent Monte Carlo treatment 
of the global dynamics and full ste llar and binary evolu- 
tion. We confirm the predictions of ISigurdsson fc Phinnevl 
1 19931 ) that BH-BH, but not NS-NS or NS-BH, binaries are 
produced efficiently in star clusters by dynamical interac- 
tions. We find that the BHs are mass-segregated and interact 
strongly with each other , confi rming the more approximate 
models of lO'Learv et all |2006l ). We find no BH-BH merg- 
ers within the clusters we simulate but many hard BH-BH 
escapers that will merge in the galactic field within a Hub- 
ble time. Detection rates for both ground- and space-based 
detectors will be calculated and presented in a future pa- 
per. We find that our simulations produce potential LISA 
sources that, while rare, will be highly eccentric and may 
still represent a significant enhancement to the galactic field 
population. We will certainly produce more detections when 
the escapers are included. We conclude that star clusters 
produce BH-BH binaries efficiently and must be taken into 
account when considering detection rates for both ground- 
and space-based detectors. 
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